/*

Preschool Availability and Female Labor Force Participation: Evidence from Indonesia
Daniel Halim, Hillary C. Johnson, Elizaveta Perova

The World Bank
East Asia Pacific Gender Innovation Lab (EAPGIL)

Apr 2020


Objective: produce tables in the paper version January 2021

*/


** graph
set scheme burd

***************************************************************************
*Table 1: DDD  validity - Test of parallel trends
***************************************************************************
use "$clean/preschool", clear
qui reg work_om nkindergov
gen used=e(sample) //
keep if year<=2002
keep if inrange(age,15,45)
global nchild "kid02 kid712 kid1318"

gen prekid=kid36_
la var prekid "has preschool kids"
gen cut=.
gen growthcut=.
gen growthcutprekid=.
gen growthprekid=.
gen cutprekid=.

**** PUBLIC
est clear
foreach v in work_om {
foreach x in gov {
foreach y in 1993 1996 1999 {
	capture ren gkinder`x' gkinder
	replace cut=year>=`y'
	replace growthcut=gkinder*cut
	replace growthcutprekid=gkinder*cut*prekid
	replace growthprekid=gkinder*prekid
	replace cutprekid=cut*prekid
	la var cut "after cut-off"
	la var growthcut "growth x after cut-off"
	la var growthprekid "growth x has preschool kids"
	la var cutprekid "after cut-off x has preschool kids"
	la var growthcutprekid "growth x after cut-off x has preschool kids"

	* 1) DDD + trend
	eststo: reg `v' growthcutprekid cutprekid growthprekid growthcut prekid $nchild i.kabcurrent i.year i.age urban  if used, cluster(kabcurrent)
	sum `v' if e(sample)
	estadd scalar varmean = r(mean)
	estadd local kabfe "X"
	estadd local yearfe "X"

	capture ren gkinder gkinder`x'
}
}
}
esttab using "$results/HJP_Table_1.rtf" ///
	, replace onecell cells(b(fmt(3) star) & se(fmt(3) par)) keep(growthcutprekid cutprekid growthprekid growthcut prekid) ///
	coeflabels(highgov "High growth public preschools") ///
	refcat(growthcutprekid "", nolabel) label nonotes compress ///
	collabels(none) modelwidth(7) ///
    mgroups("88-92 vs 93-02" "88-95 vs 96-02" "88-98 vs 99-02", pattern(1  1  1)) ///
	stats(N varmean kabfe yearfe kabtrend, fmt(%9.0fc %9.3f) ///
	labels("Observations" "Mean" "District FE" "Year FE" "District trend")) ///
	title("Table 1. Test of parallel trends prior to preschool expansion") ///
	starlevels(* 0.10 ** 0.05 *** 0.01)

***************************************************************************
*** Table 2: Effect to likelihood to have preschool aged child
***************************************************************************
use "$clean/preschool", clear
qui reg work_om nkindergov
gen used=e(sample) //
keep if inrange(age,15,45)
global nchild "kid02 kid712 kid1318"

** set panel
encode pidlink, g(numid)
xtset numid year
bys numid: gen nyear=[_N]
sum nyear

est clear
**** PUBLIC	
	
	* DDD
	eststo: qui reg kid36_ nkindergov  $nchild i.kabcurrent i.year i.age urban  if used, cluster(kabcurrent)
	sum kid36_ if e(sample)
	estadd scalar varmean = r(mean)
	estadd local kabfe "X"
	estadd local yearfe "X"

	* create lag variable
	gen Lnkindergov=L.nkindergov
	la var Lnkindergov "(Lagged) Public preschool"

	* DDD
	eststo: qui reg kid36_ Lnkindergov $nchild i.kabcurrent i.year i.age urban  if used, cluster(kabcurrent)
	sum kid36_ if e(sample)
	estadd scalar varmean = r(mean)
	estadd local kabfe "X"
	estadd local yearfe "X"


esttab using "$results/HJP_Table_2.rtf" ///
	, replace onecell cells(b(fmt(3) star) & se(fmt(3) par)) keep(nkindergov Lnkindergov) ///
	refcat(nkindergov "", nolabel) label nonotes compress ///
	collabels(none) modelwidth(7) ///
	stats(N varmean kabfe yearfe, fmt(%9.0fc %9.3f) ///
	labels("Observations" "Mean" "District FE" "Year FE")) ///
	title("Table 2. Does pre-school density affect the likelihood of having pre-school children?") ///
	starlevels(* 0.10 ** 0.05 *** 0.01)

***************************************************************************
* Table 3: Sorting
***************************************************************************
use "$clean/preschool", clear
preserve
keep if inrange(age,15,45)
keep if podes
collapse (count) all=sex (sum) eligible=kid36_ (firstnm) kinderall kinderpvt kindergov nkinderall nkinderpvt nkindergov, by(kabcurrent year)
sort kabcurrent year
by kabcurrent: egen nobs = count(year)

** arbitrary: drop kab with <= 3 obs
drop if nobs<=3

** obtain diff
foreach v in all eligible kinderall kinderpvt kindergov nkinderall nkinderpvt nkindergov {
	gen d_`v' = `v'[_n]-`v'[_n-1]
	gen d2_`v' = `v'[_n-1]-`v'[_n-2]
}

la var d_all 		"d.migration (all sample)"
la var d_eligible 	"Net migration of mothers with preschool aged kids (eligible)"
la var d_nkinderall	"Net change in all preschool density"
*la var d_nkinderpvt	"Net change in private preschool density"
la var d_nkindergov	"Net change in public preschool density"
la var d2_nkinderall "(Lagged) Net change in all preschool density"
*la var d2_nkinderpvt "(Lagged) Net change in private preschool density"
la var d2_nkindergov "(Lagged) Net change in public preschool density"


*** OLS
eststo clear

eststo: reg d_eligible d_nkindergov i.year i.kabcurrent , cluster(kabcurrent)
qui sum d_eligible if e(sample)
estadd scalar varmean = r(mean)
qui sum d_nkindergov if e(sample)
estadd scalar xvarmean = r(mean)
estadd local kabfe "X"
estadd local yearfe "X"

eststo: reg d_eligible d2_nkindergov i.year i.kabcurrent , cluster(kabcurrent)
qui sum d_eligible if e(sample)
estadd scalar varmean = r(mean)
qui sum d_nkindergov if e(sample)
estadd scalar xvarmean = r(mean)
estadd local kabfe "X"
estadd local yearfe "X"


esttab using "$results/HJP_Table_3.rtf" ///
	, replace onecell cells(b(fmt(3) star) & se(fmt(3) par)) keep(d_* d2_*) ///
	order(d_* d2_*) label nonotes compress ///
	collabels(none) modelwidth(7) ///
	title("Table 3. Do mothers sort to high preschool access districts?") ///
	stats(N varmean xvarmean kabfe yearfe, fmt(%9.0fc %9.3f) ///
	labels("Observations" "Mean" "Mean of net change in preschool density" "District FE" "Year FE")) ///
	starlevels(* 0.10 ** 0.05 *** 0.01)


***************************************************************************
* Table 4: Work status
***************************************************************************
use "$clean/preschool", clear
keep if inrange(age,15,45)
global nchild "kid02 kid712 kid1318"
*************************************** Work status

ren work_om y1
ren sidejob_om y2
ren selfemp_om y3
ren govwork_om y4
ren pvtwork_om y5
ren famwork_om y6

la var y1 "Employed"
la var y2 "Has a second job"
la var y3 "Self-employed"
la var y4 "Government worker"
la var y5 "Private worker"
la var y6 "Unpaid family worker"


***************** obtain q-value & bonferroni
*** uncorrected P-values
capt drop uP_tk_kidgov36
capt drop bonfe_tk_kidgov36
capt drop hochbergP_tk_kidgov36
gen uP_tk_kidgov36= .
gen bonfe_tk_kidgov36= .
*** corrected q-values & bonferroni
foreach j of numlist 1(1)6 {
	*** DDD
	reg y`j' i.age tk_kidgov36 nkindergov kid36_ $nchild i.year i.age urban i.kabcurrent, cluster(kabcurrent)
	replace uP_tk_kidgov36 = 2*ttail(e(N),abs(_b[tk_kidgov36]/_se[tk_kidgov36])) in `j'
	replace bonfe_tk_kidgov36=6*(2*ttail(e(N),abs(_b[tk_kidgov36]/_se[tk_kidgov36]))) in `j'
	replace bonfe_tk_kidgov36=1 if bonfe_tk_kidgov36>1 in `j'
}
qqvalue uP_tk_kidgov36, method(hochberg) qvalue(hochbergP_tk_kidgov36)


***************** RUN regression again to be reported in table
est clear
foreach j of numlist 1(1)6{
	
	eststo: reg y`j' i.age tk_kidgov36 nkindergov kid36_ $nchild i.year i.age urban i.kabcurrent, cluster(kabcurrent)
	sum y`j' if e(sample)
	estadd scalar varmean = r(mean)
	* p value
	local tstat = _b[tk_kidgov36]/_se[tk_kidgov36]
	estadd scalar pval = 2*ttail(e(df_r),abs(`tstat'))
	* fdr p value
	sum hochbergP_tk_kidgov36 in `j'
	estadd scalar fdr = r(mean)
	* bonferroni
	sum bonfe_tk_kidgov36 in `j'
	estadd scalar bonfe= r(mean)
	* fe
	estadd local kabfe "X"
	estadd local yearfe "X"

}

esttab using "$results/HJP_Table_4.rtf" ///
	, replace onecell cells(b(fmt(3) star) & se(fmt(3) par)) keep(tk_kidgov36 nkindergov kid36_) label nonotes compress ///
	collabels(none) refcat(tk_gov36 "", nolabel) modelwidth(7) ///
	title("Table 4. Effect of preschool availability on females work status") ///
	stats(fdr bonfe N varmean kabfe yearfe, fmt(%9.3f %9.3f %9.0fc %9.3f) labels("FDR q-value" "Bonferroni corrected p-values" "Observations" "Mean"  "District FE" "Year FE")) ///
	starlevels(* 0.10 ** 0.05 *** 0.01)

ren y1 work_om 
ren y2 sidejob_om 
ren y3 selfemp_om 
ren y4 govwork_om 
ren y5 pvtwork_om 
ren y6 famwork_om 

***************************************************************************
*Table 5: Earnings and hours
***************************************************************************
use "$clean/preschool", clear
keep if inrange(age,15,45)
global nchild "kid02 kid712 kid1318"
*************************************** Earning and Hours

ren lrwage_m ly1
ren lrprofit_m ly2
ren lrincome_m ly3
ren lworkhour ly4

la var ly1 "Salary"
la var ly2 "Net Profit"
la var ly3 "Income"
la var ly4 "Working hours"

la var tk_kidgov36 "Public preschool * Eligible"
la var nkinderall "Public preschools" 
la var kid36_ "Eligible child"

ren wage_m y1
ren profit_m y2
ren income_m y3
ren workhour y4


***************** obtain q-value & bonferroni
*** uncorrected P-values
capt drop uP_tk_kidgov36
capt drop hochbergP_tk_kidgov36
capt drop bonfe_tk_kidgov36
gen uP_tk_kidgov36 = .
gen bonfe_tk_kidgov36=.

*** corrected q-values & bonferroni
foreach j of numlist 1(1)4 {
	
	qui reg ly`j' i.age tk_kidgov36 nkindergov kid36_ $nchild i.year i.age urban /*if firstkid36_==1 | firstkid02_==1*/, cluster(kabcurrent)
	replace uP_tk_kidgov36 = 2*ttail(e(N),abs(_b[tk_kidgov36]/_se[tk_kidgov36])) in `j'
	replace bonfe_tk_kidgov36=4*(2*ttail(e(N),abs(_b[tk_kidgov36]/_se[tk_kidgov36]))) in `j'
	replace bonfe_tk_kidgov36=1 if bonfe_tk_kidgov36>1 in `j'
}
qqvalue uP_tk_kidgov36, method(hochberg) qvalue(hochbergP_tk_kidgov36)

***************** RUN regression again to be reported in table
est clear
foreach j of numlist 1(1)4{
	*** DDD
	eststo: qui reg ly`j' i.age tk_kidgov36 nkindergov kid36_ $nchild i.year i.age urban /*if firstkid36_==1 | firstkid02_==1*/, cluster(kabcurrent)
	sum y`j' if e(sample)
	estadd scalar varmean = r(mean)
	* p value
	local tstat = _b[tk_kidgov36]/_se[tk_kidgov36]
	estadd scalar pval = 2*ttail(e(df_r),abs(`tstat'))
	* fdr p value
	sum hochbergP_tk_kidgov36 in `j'
	estadd scalar fdr = r(mean)
	* bonferroni
	sum bonfe_tk_kidgov36 in `j'
	estadd scalar bonfe= r(mean)

	estadd local kabfe "X"
	estadd local yearfe "X"

}

esttab using "$results/HJP_Table_5.rtf" ///
	, replace onecell cells(b(fmt(3) star) & se(fmt(3) par)) keep(tk_kidgov36 nkindergov kid36_) label nonotes compress ///
	collabels(none) refcat(tk_kidgov36 "", nolabel) modelwidth(7) ///
	title("Table 5. Effect of public preschool availability on females earnings and work hours") ///
	stats(fdr bonfe N varmean kabfe yearfe, fmt(%9.3f %9.3f %9.0fc %9.3f) labels("FDR q-value" "Bonferroni corrected p-values" "Observations" "Mean"  "District FE" "Year FE")) ///
	starlevels(* 0.10 ** 0.05 *** 0.01)

ren  ly1 lrwage_m
ren  ly2 lrprofit_m
ren  ly3 lrincome_m
ren  ly4 lworkhour

ren  y1 wage_m
ren  y2 profit_m
ren  y3 income_m
ren  y4 workhour

***************************************************************************
*** Table 6: Effect to maternal employment - robustness check
***************************************************************************
use "$clean/preschool", clear
qui reg work_om nkindergov
gen used=e(sample) //
global nchild "kid02 kid712 kid1318"
keep if inrange(age,15,45) // age restriction

ren work_om y1 
ren sidejob_om y2 
ren selfemp_om y3 
ren govwork_om y4 
ren pvtwork_om y5 
ren famwork_om y6 


*** DDD-FE
*** uncorrected P-values
capt drop uP_tk_kidgov36
capt drop bonfe_tk_kidgov36
capt drop hochbergP_tk_kidgov36
gen uP_tk_kidgov36= .
gen bonfe_tk_kidgov36= .
*** corrected q-values & bonferroni
foreach v of numlist 1(1)6 {

	eststo: areg y`v' tk_kidgov36 nkindergov kid36_ $nchild i.kabcurrent i.year i.age urban, a(pidlink) cluster(kabcurrent)
	replace uP_tk_kidgov36 = 2*ttail(e(N),abs(_b[tk_kidgov36]/_se[tk_kidgov36])) in `v'
	replace bonfe_tk_kidgov36=6*(2*ttail(e(N),abs(_b[tk_kidgov36]/_se[tk_kidgov36]))) in `v'
	replace bonfe_tk_kidgov36=1 if bonfe_tk_kidgov36>1 in `v'
}
qqvalue uP_tk_kidgov36, method(hochberg) qvalue(hochbergP_tk_kidgov36)

est clear
foreach v of numlist 1(1)6 {
	la var nkindergov "Preschools density"
	la var kid36_ "Eligible child"
	la var tk_kidgov36 "Preschool * Eligible"

***Infer in between
	* DDD-FE 
	eststo: areg y`v' tk_kidgov36 nkindergov kid36_ $nchild i.kabcurrent i.year i.age urban, a(pidlink) cluster(kabcurrent)
	sum y`v' if e(sample)
    estadd scalar varmean = r(mean)
	* fdr p value
	sum hochbergP_tk_kidgov36 in `v'
	estadd scalar fdr = r(mean)
	* bonferroni
	sum bonfe_tk_kidgov36 in `v'
	estadd scalar bonfe= r(mean)

	estadd local kabfe "X"
	estadd local yearfe "X"
	estadd local indfe "X"
	estadd local spec "DDD-FE"
}
esttab using "$results/HJP_Table_6.rtf" ///
	, replace onecell cells(b(fmt(3) star) & se(fmt(3) par)) keep(tk_kidgov36 nkindergov kid36_) ///
	refcat(tk_kidgov36 "{\b Panel A}", nolabel) label nonotes compress ///
	collabels(none) modelwidth(7) ///
	title("Table 6. Robustness check: Effect of preschool on maternal employment") ///
	stats(fdr bonfe N varmean kabfe yearfe indfe spec, fmt(%9.3f %9.3f %9.0fc %9.3f) ///
	labels("FDR q-value" "Bonferroni corrected p-values" "Observations" "Mean" "District FE" "Year FE" "Individual FE" "Model")) ///
	starlevels(* 0.10 ** 0.05 *** 0.01)

***Linear projection
*** uncorrected P-values
capt drop uP_tk_kidgov36
capt drop bonfe_tk_kidgov36
capt drop hochbergP_tk_kidgov36
gen uP_tk_kidgov36= .
gen bonfe_tk_kidgov36= .
*** corrected q-values & bonferroni
foreach v of numlist 1(1)6 {

	eststo: reg y`v' tkdens_kidgov36 tkdensgov kid36_ $nchild i.kabcurrent i.year i.age urban, cluster(kabcurrent)
	replace uP_tk_kidgov36 = 2*ttail(e(N),abs(_b[tkdens_kidgov36]/_se[tkdens_kidgov36])) in `v'
	replace bonfe_tk_kidgov36=6*(2*ttail(e(N),abs(_b[tkdens_kidgov36]/_se[tkdens_kidgov36]))) in `v'
	replace bonfe_tk_kidgov36=1 if bonfe_tk_kidgov36>1 in `v'
}
qqvalue uP_tk_kidgov36, method(hochberg) qvalue(hochbergP_tk_kidgov36)

est clear
foreach v of numlist 1(1)6 {
	la var tkdensgov "Preschools density"
	la var kid36_ "Eligible child"
	la var tkdens_kidgov36 "Preschool * Eligible"

	* DDD - Linear Projection
	eststo: reg y`v' tkdens_kidgov36 tkdensgov kid36_ $nchild i.kabcurrent i.year i.age urban, cluster(kabcurrent)
	sum y`v' if e(sample)
    estadd scalar varmean = r(mean)
	* fdr p value
	sum hochbergP_tk_kidgov36 in `v'
	estadd scalar fdr = r(mean)
	* bonferroni
	sum bonfe_tk_kidgov36 in `v'
	estadd scalar bonfe= r(mean)

	estadd local kabfe "X"
	estadd local yearfe "X"
	estadd local indfe "X"
	estadd local spec "DDD- Linear projection"
}

esttab using "$results/HJP_Table_6.rtf" ///
	, append onecell cells(b(fmt(3) star) & se(fmt(3) par)) keep(tkdens_kidgov36 tkdensgov kid36_) ///
	refcat(tkdens_kidgov36 "{\b Panel B}", nolabel) label nonotes compress ///
	collabels(none) modelwidth(7) ///
	title("") ///
	stats(fdr bonfe N varmean kabfe yearfe indfe spec, fmt(%9.3f %9.3f %9.0fc %9.3f) ///
	labels("FDR q-value" "Bonferroni corrected p-values" "Observations" "Mean" "District FE" "Year FE" "Individual FE" "Model")) ///
	starlevels(* 0.10 ** 0.05 *** 0.01)

*** PODES only
*** uncorrected P-values
capt drop uP_tk_kidgov36
capt drop bonfe_tk_kidgov36
capt drop hochbergP_tk_kidgov36
gen uP_tk_kidgov36= .
gen bonfe_tk_kidgov36= .
*** corrected q-values & bonferroni
foreach v of numlist 1(1)6 {

	eststo: reg y`v' tk_kidgov36 nkindergov kid36_ $nchild i.kabcurrent i.year i.age urban if podes, cluster(kabcurrent)
	replace uP_tk_kidgov36 = 2*ttail(e(N),abs(_b[tk_kidgov36]/_se[tk_kidgov36])) in `v'
	replace bonfe_tk_kidgov36=6*(2*ttail(e(N),abs(_b[tk_kidgov36]/_se[tk_kidgov36]))) in `v'
	replace bonfe_tk_kidgov36=1 if bonfe_tk_kidgov36>1 in `v'
}
qqvalue uP_tk_kidgov36, method(hochberg) qvalue(hochbergP_tk_kidgov36)

est clear
foreach v of numlist 1(1)6 {
	la var nkindergov "Preschools density"
	la var kid36_ "Eligible child"
	la var tk_kidgov36 "Preschool * Eligible"

	* DDD PODES
	eststo: reg y`v' tk_kidgov36 nkindergov kid36_ $nchild i.kabcurrent i.year i.age urban if podes, cluster(kabcurrent)
	sum y`v' if e(sample)
    estadd scalar varmean = r(mean)
	* fdr p value
	sum hochbergP_tk_kidgov36 in `v'
	estadd scalar fdr = r(mean)
	* bonferroni
	sum bonfe_tk_kidgov36 in `v'
	estadd scalar bonfe= r(mean)

	estadd local kabfe "X"
	estadd local yearfe "X"
	estadd local indfe "X"
	estadd local spec "DDD- PODES only"
}

esttab using "$results/HJP_Table_6.rtf" ///
	, append onecell cells(b(fmt(3) star) & se(fmt(3) par)) keep(tk_kidgov36 nkindergov kid36_) ///
	refcat(tk_kidgov36 "\b Panel C", nolabel) label nonotes compress ///
	collabels(none) modelwidth(7) ///
	title("") ///
	stats(fdr bonfe N varmean kabfe yearfe indfe spec, fmt(%9.3f %9.3f %9.0fc %9.3f) ///
	labels("FDR q-value" "Bonferroni corrected p-values" "Observations" "Mean" "District FE" "Year FE" "Individual FE" "Model")) ///
	starlevels(* 0.10 ** 0.05 *** 0.01)

ren y1 work_om 
ren y2 sidejob_om 
ren y3 selfemp_om 
ren y4 govwork_om 
ren y5 pvtwork_om 
ren y6 famwork_om 

***************************************************************************
*** Table 7: Robustness check: Effect to maternal employment
***************************************************************************
use "$clean/preschool", clear
qui reg work_om nkindergov
gen used=e(sample) //
global nchild "kid02 kid712 kid1318"
keep if inrange(age,15,45) // age restriction

ren lrwage_m y1
ren lrprofit_m y2
ren lrincome_m y3
ren lworkhour y4

la var y1 "Salary"
la var y2 "Net Profit"
la var y3 "Income"
la var y4 "Working hours"

la var tk_kidgov36 "Public preschool * Eligible"
la var nkinderall "Public preschools" 
la var kid36_ "Eligible child"

ren wage_m r1
ren profit_m r2
ren income_m r3
ren workhour r4

*** DDD-FE
*** uncorrected P-values
capt drop uP_tk_kidgov36
capt drop bonfe_tk_kidgov36
capt drop hochbergP_tk_kidgov36
gen uP_tk_kidgov36= .
gen bonfe_tk_kidgov36= .
*** corrected q-values & bonferroni
foreach v of numlist 1(1)4 {

	eststo: areg y`v' tk_kidgov36 nkindergov kid36_ $nchild i.kabcurrent i.year i.age urban, a(pidlink) cluster(kabcurrent)
	replace uP_tk_kidgov36 = 2*ttail(e(N),abs(_b[tk_kidgov36]/_se[tk_kidgov36])) in `v'
	replace bonfe_tk_kidgov36=4*(2*ttail(e(N),abs(_b[tk_kidgov36]/_se[tk_kidgov36]))) in `v'
	replace bonfe_tk_kidgov36=1 if bonfe_tk_kidgov36>1 in `v'
}
qqvalue uP_tk_kidgov36, method(hochberg) qvalue(hochbergP_tk_kidgov36)

est clear
foreach v of numlist 1(1)4 {
	la var nkindergov "Preschools density"
	la var kid36_ "Eligible child"
	la var tk_kidgov36 "Preschool * Eligible"

	* DDD-FE 
	eststo: areg y`v' tk_kidgov36 nkindergov kid36_ $nchild i.kabcurrent i.year i.age urban, a(pidlink) cluster(kabcurrent)
    sum r`v' if e(sample) 
    estadd scalar varmean = r(mean)
	* fdr p value
	sum hochbergP_tk_kidgov36 in `v'
	estadd scalar fdr = r(mean)
	* bonferroni
	sum bonfe_tk_kidgov36 in `v'
	estadd scalar bonfe= r(mean)

	estadd local kabfe "X"
	estadd local yearfe "X"
	estadd local indfe "X"
	estadd local spec "DDD-FE"
}

esttab using "$results/HJP_Table_7.rtf" ///
	, replace onecell cells(b(fmt(3) star) & se(fmt(3) par)) keep(tk_kidgov36 nkindergov kid36_) ///
	refcat(tk_kidgov36 "{\b Panel A}", nolabel) label nonotes compress ///
	collabels(none) modelwidth(7) ///
	title("Table 7. Robustness check: Effect of preschool on earnings and work hours") ///
	stats(fdr bonfe N varmean kabfe yearfe indfe spec, fmt(%9.3f %9.3f %9.0fc %9.3f) ///
	labels("FDR q-value" "Bonferroni corrected p-values" "Observations" "Mean" "District FE" "Year FE" "Individual FE" "Model")) ///
	starlevels(* 0.10 ** 0.05 *** 0.01)

***Linear projection
*** uncorrected P-values
capt drop uP_tk_kidgov36
capt drop bonfe_tk_kidgov36
capt drop hochbergP_tk_kidgov36
gen uP_tk_kidgov36= .
gen bonfe_tk_kidgov36= .
*** corrected q-values & bonferroni
foreach v of numlist 1(1)4 {

	eststo: reg y`v' tkdens_kidgov36 tkdensgov kid36_ $nchild i.kabcurrent i.year i.age urban, cluster(kabcurrent)
 	replace uP_tk_kidgov36 = 2*ttail(e(N),abs(_b[tkdens_kidgov36]/_se[tkdens_kidgov36])) in `v'
	replace bonfe_tk_kidgov36=4*(2*ttail(e(N),abs(_b[tkdens_kidgov36]/_se[tkdens_kidgov36]))) in `v'
	replace bonfe_tk_kidgov36=1 if bonfe_tk_kidgov36>1 in `v'
}
qqvalue uP_tk_kidgov36, method(hochberg) qvalue(hochbergP_tk_kidgov36)

est clear
foreach v of numlist 1(1)4 {
	la var tkdensgov "Preschools density"
	la var kid36_ "Eligible child"
	la var tkdens_kidgov36 "Preschool * Eligible"

	* DDD - Linear Projection
	eststo: reg y`v' tkdens_kidgov36 tkdensgov kid36_ $nchild i.kabcurrent i.year i.age urban, cluster(kabcurrent)
    sum r`v' if e(sample)
    estadd scalar varmean = r(mean) 
	* fdr p value
	sum hochbergP_tk_kidgov36 in `v'
	estadd scalar fdr = r(mean)
	* bonferroni
	sum bonfe_tk_kidgov36 in `v'
	estadd scalar bonfe= r(mean)
	estadd local kabfe "X"
	estadd local yearfe "X"
	estadd local indfe "-"
	estadd local spec "DDD- Linear projection"
}

esttab using "$results/HJP_Table_7.rtf" ///
	, append onecell cells(b(fmt(3) star) & se(fmt(3) par)) keep(tkdens_kidgov36 tkdensgov kid36_) ///
	refcat(tkdens_kidgov36 "{\b Panel B}", nolabel) label nonotes compress ///
	collabels(none) modelwidth(7) ///
	title("") ///
	stats(fdr bonfe N varmean kabfe yearfe indfe spec, fmt(%9.3f %9.3f %9.0fc %9.3f) ///
	labels("FDR q-value" "Bonferroni corrected p-values" "Observations" "Mean" "District FE" "Year FE" "Individual FE" "Model")) ///
	starlevels(* 0.10 ** 0.05 *** 0.01)

*** PODES only
*** uncorrected P-values
capt drop uP_tk_kidgov36
capt drop bonfe_tk_kidgov36
capt drop hochbergP_tk_kidgov36
gen uP_tk_kidgov36= .
gen bonfe_tk_kidgov36= .
*** corrected q-values & bonferroni
foreach v of numlist 1(1)4 {

	eststo: reg y`v' tk_kidgov36 nkindergov kid36_ $nchild i.kabcurrent i.year i.age urban if podes, cluster(kabcurrent)
	replace uP_tk_kidgov36 = 2*ttail(e(N),abs(_b[tk_kidgov36]/_se[tk_kidgov36])) in `v'
	replace bonfe_tk_kidgov36=4*(2*ttail(e(N),abs(_b[tk_kidgov36]/_se[tk_kidgov36]))) in `v'
	replace bonfe_tk_kidgov36=1 if bonfe_tk_kidgov36>1 in `v'
}
qqvalue uP_tk_kidgov36, method(hochberg) qvalue(hochbergP_tk_kidgov36)

est clear
foreach v of numlist 1(1)4 {
	la var nkindergov "Preschools density"
	la var kid36_ "Eligible child"
	la var tk_kidgov36 "Preschool * Eligible"

	* DDD PODES
	eststo: reg y`v' tk_kidgov36 nkindergov kid36_ $nchild i.kabcurrent i.year i.age urban if podes, cluster(kabcurrent)
    sum r`v' if e(sample) 
    estadd scalar varmean = r(mean)
	* fdr p value
	sum hochbergP_tk_kidgov36 in `v'
	estadd scalar fdr = r(mean)
	* bonferroni
	sum bonfe_tk_kidgov36 in `v'

	estadd scalar bonfe= r(mean)
	estadd local kabfe "X"
	estadd local yearfe "X"
	estadd local indfe "-"
	estadd local spec "DDD- PODES only"
}

esttab using "$results/HJP_Table_7.rtf" ///
	, append onecell cells(b(fmt(3) star) & se(fmt(3) par)) keep(tk_kidgov36 nkindergov kid36_) ///
	refcat(tk_kidgov36 "\b Panel C", nolabel) label nonotes compress ///
	collabels(none) modelwidth(7) ///
	title("") ///
	stats(fdr bonfe N varmean kabfe yearfe indfe spec, fmt(%9.3f %9.3f %9.0fc %9.3f) ///
	labels("FDR q-value" "Bonferroni corrected p-values" "Observations" "Mean" "District FE" "Year FE" "Individual FE" "Model")) ///
	starlevels(* 0.10 ** 0.05 *** 0.01)

ren  y1 lrwage_m
ren  y2 lrprofit_m
ren  y3 lrincome_m
ren  y4 lworkhour

ren  r1 wage_m
ren  r2 profit_m
ren  r3 income_m
ren  r4 workhour

***************************************************************************
*** Appendix Table 1: Summary statistics
***************************************************************************
use "$clean/preschool", clear
keep if inrange(age,15,45)
*************************** Individual

*** label variables
la var age "Age"
la var kid36_ "Have preschool-aged child"
la var work_om "Work participation"
la var nkindergov "Public preschool density (Inferred in-between)"
la var nkinderpvt "Private preschool density (Inferred in-between)"
la var tkdensgov "Public preschool density (Linear projection)"
la var tkdenspvt "Private preschool density (Linear projection)"
la var urban "Urban"

recode urban (2=0) // 2=rural

***** store summary stats
est clear
* multiple obs per individual
estpost su age kid36_ work_om nkindergov nkinderpvt tkdensgov tkdenspvt cdens urban

***** table: all stacked
esttab using "$results/HJP_Table_1_Ind_Multiple.rtf" ///
	, replace cells("count(fmt(%12.0fc)) mean(fmt(%12.2fc)) sd(fmt(%12.2fc))") ///
	collabels(Obs Mean SD) label  nonumber noobs

*************************** Individual (PODES only)
***** store summary stats
est clear
* multiple obs per individual
estpost su age kid36_ work_om nkindergov nkinderpvt cdens urban if podes

***** table: all stacked
esttab using "$results/HJP_AppTable_1_Ind_Multiple_PODES.rtf" ///
	, replace cells("count(fmt(%12.0fc)) mean(fmt(%12.2fc)) sd(fmt(%12.2fc))") ///
	collabels(Obs Mean SD) label  nonumber noobs


**************** 1 obs per individual
egen tagid = tag(pidlink)
unique pidlink // # mothers = 10340

merge m:1 pidlink using "$ifls_clean1/b4_kw2", keep(1 3) ///
	keepusing(firstwed_yr firstwed_age)
drop _merge
merge m:1 pidlink using "$ifls_clean2/b4_kw2", update keep(1 3 4 5) ///
	keepusing(firstwed_yr firstwed_age)
drop _merge
merge m:1 pidlink using "$ifls_clean3/b4_kw2", update keep(1 3 4 5) ///
	keepusing(firstwed_yr firstwed_age)
drop _merge
merge m:1 pidlink using "$ifls_clean4/b4_kw2", update keep(1 3 4 5) ///
	keepusing(firstwed_yr firstwed_age)
drop _merge
merge m:1 pidlink using "$ifls_clean5/b4_kw2", update keep(1 3 4 5) ///
	keepusing(firstwed_yr firstwed_age)
drop _merge

* define age of first marriage
gen marryage = firstwed_age
replace marryage = firstwed_yr - bth_year if marryage==.
drop firstwed_yr firstwed_age
gen marry_year = bth_year + marryage

*** generate variables
egen obspermother = count(year), by(pidlink)
egen podespermother = sum(podes), by(pidlink)

*** label variables
la var numivw "Number of surveys"
la var obspermother "Number of years"
la var podespermother "Number of PODES years"
la var age_birth1 "Age of first birth"
la var marryage "Age of first marriage"
la var child "Number of children"
la var modeyos "Years of education"

***** store summary stats
est clear
* 11 obs per individual
estpost su numivw obspermother podespermother marryage age_birth1 child modeyos if tagid==1

***** table: all stacked
esttab using "$results/HJP_AppTable_1_Ind_Unique.rtf" ///
	, replace cells("count(fmt(%12.0fc)) mean(fmt(%12.2fc)) sd(fmt(%12.2fc))") ///
	collabels(Obs Mean SD) label  nonumber noobs



*************************** District
use "$clean/PODES/kindergarten across podes",clear

*** reshape
ren kabcode_des93 kab
ren *90 *1990
ren *93 *1993
ren *96 *1996
ren *00 *2000
ren *03 *2003
ren *05 *2005
ren *08 *2008
ren *11 *2011
ren *14 *2014

keep kindergov* kinderpvt* cdens* kab
reshape long kindergov kinderpvt cdens, i(kab) j(year)

*** define kindergarten availability per 1000 children
foreach v in gov pvt {
	* define density
	gen dens`v' = kinder`v'/cdens*1000
}

*** label variables
la var kindergov "Public preschool count"
la var kinderpvt "Private preschool count"
la var densgov "Public preschool density"
la var denspvt "Private preschool density"
la var cdens "Child population"

***** store summary stats
est clear
estpost su kindergov kinderpvt cdens densgov denspvt

***** table: all stacked
esttab using "$results/HJP_AppTable_1_Kab.rtf" ///
	, replace cells("count(fmt(%12.0fc)) mean(fmt(%12.2fc)) sd(fmt(%12.2fc))") ///
	refcat(kindergov "{\b Number of Districts = 290}", nolabel) ///
	collabels(Obs Mean SD) compress label  nonumber noobs